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We study the structure of neutron stars in f(R) gravity theories with perturbative constraints. 
We derive the modified Tolman-Oppenheimer-Volkov equations and solve them for a polytropic 
equation of state. We investigate the resulting modifications to the masses and radii of neutron 
stars and show that observations of surface phenomena alone cannot break the degeneracy between 
altering the theory of gravity versus choosing a different equation of state of neutron-star matter. 
On the other hand, observations of neutron-star cooling, which depends on the density of matter at 
the stellar interior, can place significant constraints on the parameters of the theory. 



I. INTRODUCTION 

Recent interest in modified theories of gravity has been 
spurred by the discovery that the Universe is undergo- 
ing accelerated expansion (see, e.g., [IHl])- The simplest 
solution consistent with these observations posits a cos- 
mological constant A. The magnitude of this cosmologi- 
cal constant is significantly less than what was expected, 
and many undertakings have been made to see if there 
are plausible alternative explanations [J, Q. Outstand- 
ing questions also present themselves in the formation of 
singularities @ and the seeming contradiction between 
quantum mechanics and gravity in the context of black 
hole thermodynamics Q. All these suggest that there 
may yet be much to understand about the nature of grav- 
ity at extreme-curvature scales, far removed from our ev- 
eryday experience. 

The two most popular approaches to modifying grav- 
ity have been the introduction of an additional scalar 
field (e.g. 0), or the related approach of replacing the 
Einstcin-Hilbert action with a general function of the 
Ricci scalar f(R) (e.g. Q). Within either framework the 
additional scalar degree of freedom can be tuned to mimic 
the cosmological constant, or any type of cosmological 
evolution at cosmological scales @ . 

Despite the premise of such modifications, the non- 
linear character of gravitational theories has proven a 
significant obstacle to introducing new dynamical fields 
to drive modifications to gravity at the cosmological scale 
without the same fields reemerging at widely different 
curvature scales. One such example is the problem of 
ensuring that f(R) = i?± /z 4 /R theories pass the current 
Parametrized Post-Newtonian (PPN) bounds. When the 
new field is dynamical, the PPN parameter 7 is forced to 
a value of 1/2, which is very far from the present ex- 
perimental bound [lOj . As a result one has to choose a 
function f(R) only from the class which can adequately 
suppress the new dynamical field on solar-system scales. 



The chameleon mechanism [llMl3l | provides such an al- 
ternative. 

In addition to the PPN constraints, instabilities related 
to the functional form of f(R) have also been studied at 
length. This is especially true for the Dolgov-Kawasaki 
instability [l4j], which requires that d 2 f/dR 2 > in order 
that the effective mass of the equivalent scalar degree of 
freedom be positive. In the strong-field regime, recent re- 
sults [15[ suggest that this very choice may well prohibit 
the formation of compact objects above a curvature scale 
readily observed. However, the fatal curvature singular- 
ity may be avoided by the chameleon mechanism [li* [l7| . 

Perhaps the source of the instabilities and consistency 
issues many of these models encounter is the result of 
treating these modifications as though they are exact. 
The original motivation behind introducing additional 
functions of the curvature was to generate a new phe- 
nomenology at a specific scale. However, many of the 
problems encountered by f(R) gravity theories originate 
at curvature scales far removed from the ones under con- 
sideration. An alternative formulation for handling cor- 
rections to General Relativity is to view the new terms 
as only the next to leading order terms in a larger expan- 
sion. In this context there is no reason to suspect that 
the new phenomenology is due to new dynamical fields. 
The technique for handling a field expansion of this form 
is well developed [l8| and is known as perturbative con- 
straints or order reduction (l9j . 

Gravity with perturbative constraints allows us to ex- 
plore alternative phenomenologies of gravity while main- 
taining important consistency conditions including gauge 
invariance, the assumption that we are approximating a 
fundamentally second order field theory, and the con- 
servation of stress-energy. Maintaining such constraints 
while enlarging the space of possible behaviors of gravi- 
tation is t he g oal also of the Parametrized Post-Friedman 
approach |2(H22j. 

In previous works [H, [24| , we have analyzed the effect 
of treating f(R) models of gravity via perturbative con- 
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straints primarily at cosmological scales. In this paper, 
we examine the ramifications of modifications to grav- 
ity in the context of compact objects. We show how the 
method of perturbative constraints allows for a consistent 
phenomenology for gravity on both large (Hubble-length 
perturbations linear in metric variables, but strongly rel- 
ativistic, L ~ c/Hq) and small scales (stellar scales, non- 
linear in metric perturbations, and strongly relativistic, 
GM/rc 2 ~ 1.) 

The layout of this work is as follows. In Section UU we 
review the equations of f(R) gravity treated with pertur- 
bative constraints. In Section UTTl we derive the modified 
Tolman-Oppenheimer-Volkov equations and show that 
the exterior solution is the Schwartzchild-de Sitter met- 
ric. In Section HVl we demonstrate that such objects are 
stable and we solve numerically for their mass-radius re- 
lation for a polytropic equations of state. Finally in Sec- 
tion El we discuss how we can discriminate modifications 
to gravity from uncertainty in the neutron star equation 
of state. 



II. PERTURBATIVE CONSTRAINTS 



The perturbative consistency of this approach is guar- 
anteed to order n provided a n+1 g^Il +1 ^ <C g$ + . . . + 
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a" as we outlined in a previous paper |24|. Note 

that this condition is not to be understood as requiring 
the product af(R) to be necessarily smaller in magnitude 
than R. 

For the purposes of this work it will prove useful to 
rewrite Eq. @ using its trace 



R-a [f R R -2f + 3Df R ] + 0(a 2 



-8ttT + 4A . (4) 



Substituting the Ricci scalar R from the above equation 
into Eq. ([2]) gives 
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0{a 2 ) =8n[T t 



.(5) 



This is the form of the field equation we will be using. 
Henceforth we shall understand the equality sign to mean 
equality up to order a and drop the explicit use of 0(a 2 ). 



Gravity with perturbative constraints [1 81 ] (or order- 
reduction [19]) is a technique for treating equations of 
motion that appear higher than second order, where the 
origin of the higher derivatives can be traced to the trun- 
cation of an infinite series expansion. Such a situation 
can arise with non-local theories as well as effective field 
theories. 

In the context of f(R) gravity theories, we parametrize 
the deviation from General Relativity by a single param- 
eter a and derive the equation of motion from a covariant 
action 



■5' 



1 

167r 



d 4 x y/=g~ [R-2A + af(R) + 0(a 2 )] 

+S M {g^,^) , (1) 



with G = c = 1. Here g^ v is the metric, g its determi- 
nant, and R the Ricci scalar. We denote any additional 
terms above order a by 0{a 2 ). We may not impose any 
constraints at the level of the action without altering the 
nature of the variational principle. The resulting field 
equation is 



R 



1 

-g^vR + .g^A + a 



SrRhu — -^g^uf— 



(V^V, - g^U) f R + 0{a 2 ) = 8ttT^ , 



(2) 



where f R = df/dR. 

At zeroth order in a, these equations are second order 
in the metric; we denote the solution at this order by g$ ■ 
We then solve the system for the higher order terms by 
writing 



g(P) 



0{a 2 



(3) 



III. STARS WITH PERTURBATIVE 
CONSTRAINTS 

The metric of a static, spherically symmetric object 
can always be written in the form 

ds 2 = -B(r)dt 2 + A(r)dr 2 + r 2 (d6 2 + sin 2 6dcj> 2 ) , (6) 

where B{r) = (r) + aB^ (r) + . . ., A(r) = #(r) + 
aA^\r) + . . ., and B (0) (r) and yl (0) (r) are the general 
relativistic metric elements. 

For the purpose of this paper we presume the form 
f(R) oc R n+1 for an integer n ^ 0,-1. We shall also as- 
sume that the energy-momentum tensor within the star is 
that of a perfect fluid. Following our previous studies [24| 
we find it convenient to express the 0(a) correction in 
terms of the derivative f R . The first three field equations 
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where the prime denotes differentiation with respect to r. 
The fourth field equation is identical to Eq. ^ because 
of the symmetry of the spacetime. Terms with a factor 
fa preceding them are already first order in the small 
parameter a so all such terms should be evaluated at 
order 0(a ), where for example 



i? (0) = 8tt - 3P) + 4A 



and 



M (0) = 4tt 



pr' 



dr 



(10) 
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In order to motivate the form of the metric element A(r) 
that we will be using, we first examine the solution exte- 
rior to the star. 



A. The Exterior Metric 

To solve for the exterior solution to Eq. ([5]) , we require 
that outside the star T M „ = 0. Therefore, at O(a ), the 
exterior metric satisfies 



(12) 



where R$ is the Ricci tensor derived from the metric to 
O(a ). Consequently the Ricci scalar at 0(oP) is R^ = 
4A. 

Note from equations and © that the 0(a) 

correction is multiplied by a term /r oc \(n + 1)R^~\ . 
For n > 1 such a theory will allow a solution with a 
Minkowski exterior as well as solutions with A ^ 0, while 
for n < —2 the appearance of R^ in the denominator 
requires that only solutions with A ^ exist. 

In order to calculate the corrections to the vacuum so- 
lution at successively increasing orders in a, we first in- 
vestigate the perturbative term in the field equation ([5]) , 
when the Ricci curvature is constant. At 0(a) the cor- 
rection term is proportional to 

fP*# - \$ {f^R {0) f (0) ) « (» - 1) < 3 , (i3) 

where we evaluated everything explicitly in terms of R^ 
and Rpv ■ This last relation shows that, in f(R) theories 
with n = 1, the correction term in the field equation 
vanishes and hence the exterior solution is identical to 

gr m. 

We can proceed in the same manner to arbitrary orders 
in 0(a m ). The result can be formally written as 
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where the precise form of the function T is determined 
by the choice of the function f(R). 

The vacuum equations, therefore, choose a unique so- 
lution, the Schwartzchild-de Sitter metric, with 



2M Ar 2 
A(r)=[l--- — 



(15) 



and A(r)B(r) = 1 [26]. The only difference from the 
general relativistic exterior metric will be in the value of 
the effective cosmological constant, which in the case of 
f{R) gravity is 



A = J" (a, a 2 , . . . ,a m ) A 



(16) 



As a result, the PPN parameters [27| for an arbitrary 
choice of f(R) will be practically those of General Rela- 
tivity (see also discussion in [25j). 



B. Interior Solution 

In the following, we shall suppress the explicit appear- 
ance of A in the field equations by the useful redefinitions 
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(17) 



Subject to these normalizations and given the form of the 
exterior solution we shall define 



Mr) = 



1 - 



2M(r) 



(18) 



We will use this definition to all orders in the small pa- 
rameter a, with the term M{r) acquiring corrections at 
each successive order, as it is shorthand for a metric ele- 
ment. 

From the form of the elements of the Ricci tensor, and 
the above definition we obtain 

^oo fti R22 _ 2M 
~2B YA ~ ~ ~ 

Combining this with equations (JT])-© and evaluating the 
result to order 0(a) we derive the equation for mass con- 
servation in f(R) gravity with perturbative constraints 
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and arrive at the equation of hydrostatic equilibrium via 
equation ([9]) 
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Note that in solving Eqs. (|20)) and (|23j) in practice the 
evolution of the density and pressure are determined in 
terms of the familiar Tolman-Oppenheimer-Volkov equa- 
tions 



dr 



and 



dPo 
dr 



(po + Po) (M (0) + 47rP r : 



(24) 



(25) 



Here po and Po are understood to be the pressure and 
density evaluated at 0(a°), whereas we will denote the 
pressure evolved via equation (j2"3")l to 0(a) by Pi. 



IV. NUMERICAL MODELS OF NEUTRON 
STARS 

The equations we have derived so far are general and 
accommodate any choice for the correction f(R) to the 
Einstein-Hilbert action. However, in constructing numer- 
ical models of neutron stars in f(R) theories, we need to 
specify at this point the particular value of the parameter 
n we will use. 

In order to address concerns for the structure and sta- 
bility of neutron stars in cosmologically motivated mod- 
ifications of gravity (see [lj| and [Hj]), we might con- 
sider the case n = — 2 (i.e., f(R) = R^ 1 ). Since the 
matter density and pressure directly determine the Ricci 
scalar, we would anticipate such a term to be the leading 
order correction for small-curvature scales. Unlike the- 
ories with additional degrees of freedom, however, and 
as we would expect given the magnitude of R, the low- 
curvature corrections lead to no observable differences 
in the structure of compact objects. Our analysis of 
stars with these low-curvature corrections demand that 
the perturbative parameter a not be significantly larger 
than A. Such a small correction leads to no discernible 
distinction from the predictions of general relativity. 

For this reason, we will study below the case with 
n = 1, i.e., gravity theories with f(R) = R 2 . This rep- 
resents the next to leading order correction in a high- 
curvature expansion of the action. It is this regime where 
we expect the correction to be most noticeable in the case 
of compact objects. 

We choose the polytropic equation of state 
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r - i 



(26) 



for the interior of the neutron star, where T is the poly- 
tropic index. Realistic neutron star equations of state 
can be parameterized by piecewise polytropic equations 
of state [29|, H3] with r ~ 1 — 3. The lower the polytropic 
index, the stiffer the associated mass-radius relationship. 
For this study we have chosen T = 9/5 which is consistent 
with the constraints on T2 in Ref. 301. 
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FIG. 1: The mass of a neutron star as a function of its central 
density in an f(R) = R 2 gravity theory for different values of 
the small parameter a. The index of the polytropic equation 
of state was set to F = 9/5. 



We utilize the same dimensionless variables as in 
Ref. [3l|, namely 
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Because of this normalization of the various physical 
quantities, our results are independent of the normal- 
ization K of the polytropic equation of state. 

We use a fourth order Runge-Kutte integrator with 
adaptive stepsize to solve for the mass M and radius 
R of the star. We start at the center of the star by 
specifying its density (and corresponding pressure) there 
and integrate out to its surface defined where the pressure 
vanishes. 

Figure [1] shows the dependence of the mass of a neutron 
star on its central density in an f(R) — R 2 theory, for 
different values of the small parameter a. The central 
line corresponds to neutron stars in general relativity. As 
expected, for stable neutron stars the deviation from the 
general relativistic case becomes significant as the central 
density of the neutron star increases, since it is the matter 
density that directly determines the value of the Ricci 
scalar curvature. Moreover, the sign of the deviation is 
determined by the sign of the perturbative parameter 
a. By properly choosing the sign and magnitude of this 
parameter, we can cause an increase or a decrease in the 
maximum mass of stable neutron stars for a particular 
central density. We can also support stars of a certain 
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FIG. 2: The ratio £ (see Eq. I32p as a function of the param- 
eter a, for stars with central densities logp c = —2, —3, —4 
(dotted, solid and dashed lines respectively). The ratio £ mea- 
sures the degree of perturbative validity of the stellar models. 
A necessary condition for perturbative validity is £ < 1. 



mass and radius for a range of central densities and a. 

The maximum allowed magnitude of the deviations 
from the general relativistic predictions is, of course, 
constrained by the requirement that the solutions retain 
their perturbative validity. Though this constraint does 
not have a ready analytic expression, we can neverthe- 
less explore after the fact the perturbative validity of each 
stellar model. 

In particular, as a measure of the deviations from the 
general relativistic solution we choose the ratio 
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(32) 



This ratio varies with radius inside the neutron star. It 
achieves, however, its highest value at or near the center 
of the star, where the density (and hence the curvature) 
is large. Because we require the entire solution to be 
perturbatively close to the general relativistic one, we 
will evaluate the ratio £ at its maximum. A necessary 
condition for perturbative validity is £ < 1. 

Figure [2] shows the maximum ratio £ as a function of 
the parameter a. This figure demonstrates that neutron 
stars in /(i?) = R 2 theories can certainly be treated per- 
turbatively as long as —0.015 < a < 0.01. 

Of particular interest from an observational point of 
view is the mass-radius relation for neutron stars. We 
show this relation, for the same polytropic equation of 
state, in Figure[3J Depending on the value and sign of the 
parameter a, we obtain stars with larger or smaller radii 
compared to their general relativistic counterparts of the 
same gravitational mass. The extent of this variation is 
constrained by perturbative validity, which prevents the 
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FIG. 3: The mass-radius relation of neutron stars for a poly- 
tropic equation of state with F = 9/5, in an f(R) = aR 2 
gravity for different values of the parameter a. 



onset of dynamical features such as spontaneous scalar- 
ization 1321. 



V. DISCUSSION 

The predicted mass-radius relation for neutron stars 
in f(R) gravity shown above differs from that computed 
within general relativity. However, very similar devia- 
tions in the mass-radius relation can also be obtained 
within general relativity by simply changing the poly- 
tropic index of the equation of state (see [301 ] for ex- 
amples). Because the equation of state of neutron-star 
matter is weakly constrained by current experiments, 
neutron-star observables that depend only on the mass 
and radius of the star cannot distinguish between small 
differences in the equation of state versus small modifi- 
cations to gravity. 

In [33j | it was shown that observables that depend also 
on the effective surface gravity of neutron stars can break, 
in principle, this degeneracy. In particular it was shown 
that the Eddington luminosity L|? of a bursting neutron 
star depends directly on its effective surface gravity as 
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In this equation, m p is the mass of the proton, X is the 
hydrogen mass fraction in the neutron-star atmosphere, 
ctt is the Thomson scattering cross section, and 



z s = 1 - 



2M\ 1 



(34) 



is the gravitational redshift from the neutron star sur- 
face. The parameter rj is the ratio of the effective surface 
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FIG. 4: The central density p of neutron stars with different 
masses M as a function of the parameter a. Larger positive 
deviations from general relativity require larger central den- 
sities for the same neutron-star mass and, therefore, lead to 
shorter cooling times. On the other hand, larger negative de- 
viations require smaller central densities and lead to longer 
cooling time. 



gravity of the neutron star to that calculated in GR, i.e., 
,/ - (35) 
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We can calculate easily the value of the parameter r\ 
for the f{R) = R 2 theory considered here. From the 
conservation equation (|2ip we can write 
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We can then evaluate the hydrostatic equilibrium equa- 
tion (|23[) to first order in a by noting that 
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As a result equation (|38| becomes 
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At the surface layer of the neutron star p = P = and 
hence 



(40) 
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Which has the same dependence on mass and radius as 
gGR does. As a result measuring rj alone will not suffice 
to break the degeneracy due to the equation of state. 

Nevertheless constraining observationally the cooling 
rates of neutron stars can offer a discriminant. A neutron 
star cools both through photon and neutrino emission. 
The photon luminosity is determined by the temperature 
at the photosphere, which in turn depends on the den- 
sity of the photosphere. However neutrino cooling, which 
depends more sensitively on temperature than photon 
cooling does, becomes dominant for neutron stars with 
temperatures above 10 10 -ftT, and indeed is the primary 
mechanism of cooling for young neutron stars (see [34| 
for a detailed review). The high temperature and low 
interaction rate make neutrino cooling particularly sen- 
sitive to the central density of the neutron star. Figure 0] 
shows the relation between the parameter a and the cen- 
tral density of a neutron star, for three different values 
of the mass M = 0.15, 0.125, and 0.1. Large positive 
deviations from general relativity, as measured by the 
parameter a require larger central densities for neutron 
stars of a given mass, whereas the opposite is true for 
large negative deviations. As a result, because the cool- 
ing timescale scales with central density, observations of 
the surface temperatures of young neutron star can lead 
to useful constraints on the deviations from general rela- 
tivity in an f{R) gravity model, especially if the neutron- 
star masses are known. 

We will study the constraints imposed on f(R) gravity 
by current measurements of cooling rates of neutron stars 
in our galaxy in a forthcoming paper. 
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